In [1]:
from pulp import *
import pandas as pd
import matplotlib.pyplot as plt 

import warnings
warnings.filterwarnings("ignore")

Optimizando la generación de electricidad¶

No description has been provided for this image

Las empresas de gestión de la energía necesitan anticiparse a la demanda de electricidad que va a haber en un determinado período. Una vez estimada esa demanda, y teniendo en cuenta qué centrales poseen, cuál es la capacidad de generación de cada una, y qué costes de operación tienen, necesitan decidir cuánta electricidad producen con cada una de sus centrales en cada momento para garantizar que la demanda es satisfecha, a la vez que incurren en el menor coste posible.

En esta práctica vamos a simular que somos un operador de una de estas empresas, y tenemos que decidir qué hacer para gestionar la demanda de electricidad de un lugar ficticio.

Datos¶

Estimación de demanda¶

Partimos de un mes de demanda total estimada (Julio de 2023), y consideramos que vamos a optimizar justo antes de que se inicie el mes. Si leemos las demandas predichas obtenemos lo siguiente:

In [2]:
import pandas as pd
hourly_demand = pd.read_csv('./demands.csv', parse_dates=['hora'], index_col='hora').sort_index()
hourly_demand
Out[2]:
demanda
hora
2023-07-01 01:00:00 22424.578280
2023-07-01 02:00:00 20288.860640
2023-07-01 03:00:00 17962.060638
2023-07-01 04:00:00 16585.284168
2023-07-01 05:00:00 15976.342990
... ...
2023-07-31 20:00:00 27991.637100
2023-07-31 21:00:00 27942.225340
2023-07-31 22:00:00 27630.460640
2023-07-31 23:00:00 27299.872400
2023-08-01 00:00:00 26400.119460

744 rows × 1 columns

Donde vemos que tenemos una estimación por hora (31 días · 24 horas = 744 filas), medida en MW. Hay que entender que la marca de hora es la del final del periodo, es decir, que entre las 00 horas y la 01 de la noche del 1 de julio se requieren 22414 MW, mientras que entre las 23 horas del 31 de julio y la medianoche del 1 de agosto necesitamos proveer 26400 MW.

Si pintamos todo el mes, podemos comprobar que la demanda sigue un patrón de picos y valles, de forma que hay un pico por cada día (si bien la magnitud de cada pico y valle depende de factores como si es un día laborable o no, la temperatura que hace, etc.):

In [3]:
hourly_demand.plot(figsize=(16, 4), legend=False, ylabel='demanda (MW)', xlabel='día y hora');
No description has been provided for this image

Dentro de un día, también es bien conocido que hay horas valle y horas pico. En concreto, si trabajamos con el primer día de todos (1 de julio):

In [4]:
daily_demand = hourly_demand.iloc[:24]   # los 24 primeros registros
daily_demand.T   # trasponer para que no ocupe tanto en pantalla
Out[4]:
hora 2023-07-01 01:00:00 2023-07-01 02:00:00 2023-07-01 03:00:00 2023-07-01 04:00:00 2023-07-01 05:00:00 2023-07-01 06:00:00 2023-07-01 07:00:00 2023-07-01 08:00:00 2023-07-01 09:00:00 2023-07-01 10:00:00 ... 2023-07-01 15:00:00 2023-07-01 16:00:00 2023-07-01 17:00:00 2023-07-01 18:00:00 2023-07-01 19:00:00 2023-07-01 20:00:00 2023-07-01 21:00:00 2023-07-01 22:00:00 2023-07-01 23:00:00 2023-07-02 00:00:00
demanda 22424.57828 20288.86064 17962.060638 16585.284168 15976.34299 15686.931226 14952.81358 14618.695932 14675.16652 15365.754756 ... 22663.40182 23579.55476 23844.51946 24329.69594 24988.99004 25165.75476 24815.16652 24461.6371 23335.93122 22047.99004

1 rows × 24 columns

Vemos claramente que durante la madrugada la demanda baja, para después ir subiendo a lo largo de la mañana, alcanzar la demanda máxima a media tarde, y luego volver a ir decayendo según avanza la noche:

In [5]:
daily_demand.plot(marker='o', legend=False, figsize=(8, 4), ylabel='demanda (MW)');
No description has been provided for this image

Centrales¶

Lo normal es que una compañía disponga de varios tipos de centrales, cada una de ellas con unas características propias. Habitualmente la primera distinción que se hace es por el tipo de energía de cada una.

En este caso vamos a suponer que operamos 10 centrales diferentes, pertenecientes a 6 tipos de energía (carbón, nuclear, petróleo, gas, hidráulica y biomasa) y denominadas con las letras de la A la J.

A continuación figuran la fuente de energía de cada una, el combustible que consumen (que no es exactamente lo mismo que la fuente de su energía, pero que es el factor que determina los costes, como veramos más adelante), así como sus capacidades máximas de producción:

In [6]:
plants = pd.read_csv('plants.csv', index_col='central').sort_index()
plants
Out[6]:
tipo combustible producción
central
A carbón carbón sub-bituminoso 5720.22000
B carbón carbón bituminoso 5244.74856
C carbón carbón bituminoso 861.36600
D nuclear nuclear 3252.82456
E nuclear nuclear 4177.39200
F petróleo licor negro 151.67496
G biomasa licor negro 941.91944
H gas gas natural 3522.61048
I hidráulica agua 536.11236
J carbón carbón bituminoso 8793.66424

Costes¶

Dentro de los costes, a cada central le son aplicables los siguientes:

  • El del combustible que usa (es decir, lo que cuesta producir con ella 1 MWh).
  • El de operación (lo que cuesta mantenerla operativa durante una hora, independientemente de lo que esté produciendo).
  • De encendido y apagado (lo que cuesta encenderla o apagarla, que asumimos que es idéntico en ambos casos).
  • Sanitario (una estimación del coste en salud que supone usar una central; esto es solamente aplicable a las centrales de carbón A, B y C, y es proporcional a la cantidad de carbón quemado).

Los tres primeros costes se pueden considerar constantes a lo largo del mes y están recogidos en la siguiente tabla:

In [7]:
plant_costs = pd.read_csv('./plant_costs.csv', index_col='combustible').sort_index()
plant_costs
Out[7]:
coste_comb coste_oper coste_on_off
combustible
agua 0.000 8.176099 0.0
carbón bituminoso 18.050 7.384495 147.0
carbón sub-bituminoso 9.195 7.380880 147.0
gas natural 17.435 5.287164 72.0
licor negro 7.440 5.295063 72.0
nuclear 3.615 16.310160 0.0

En cambio, los costes sanitarios están calculados por cada hora individualmente. Puedes consultarlos en esta otra tabla:

In [8]:
health_costs = pd.read_csv('./health_costs.csv', parse_dates=['hora'], index_col=['hora', 'central']).sort_index()
health_costs
Out[8]:
coste_salud
hora central
2023-07-01 01:00:00 A 2.148816
B 9.030773
C 3.422946
2023-07-01 02:00:00 A 2.022208
B 7.706546
... ... ...
2023-07-31 23:00:00 B 13.274831
C 14.214563
2023-08-01 00:00:00 A 2.562029
B 12.606380
C 10.906997

2232 rows × 1 columns

Si pintamos estos costes a lo largo de todo el mes, vemos que por lo general las centrales B y C tienen mayor coste que la A, pero no forman un patrón tan claro como el de la demanda:

In [9]:
health_costs.reset_index().pivot(index='hora', columns='central', values='coste_salud').plot(figsize=(16, 4), ylabel='coste sanitario', xlabel='día y hora');
No description has been provided for this image

Restringiéndolo al día de interés (1 de julio), los costes son bastante bajos (menos de 20 por hora), pero claramente superiores en las horas nocturnas que en las diurnas:

In [10]:
daily_health_costs = health_costs.iloc[:72]   # 3 centrales, cada una de ellas 24 horas -> 72 filas
daily_health_costs.reset_index().pivot(index='hora', columns='central', values='coste_salud').plot(marker='o', figsize=(8, 4), ylabel='coste sanitario');
No description has been provided for this image

Restricciones y suposiciones¶

Existen una serie de limitaciones establecidas por razones de seguridad y eficiencia:

  • Las centrales nucleares deben estar siempre operativas (sólo se pueden apagar por lo que se conoce como "paradas programadas", cosa que ignoramos aquí).
  • Si una central está operativa (en marcha), tiene que producir al menos un cierto porcentaje de su capacidad total. Este mínimo es de:
    • Un 80% para las centrales nucleares.
    • Un 1% para el resto.
  • Además, para evitar cambios bruscos, se establece también un porcentaje máximo de cambio entre la producción de una hora y de la siguiente (ya sea por aumento o por disminución de producción de energía). En concreto, este límite es de:
    • Un 20% para las centrales nucleares.
    • Un 25% para las de carbón.
    • El 100% para el resto (es decir, que pueden pasar de estar apagadas a producir al 100% de su capacidad, o viceversa, pasar de producir al 100% a estar apagadas en la hora siguiente).

Asímismo, podemos por simplicidad asumir una serie de hechos:

  • El encendido y apagado de una central es instantáneo una vez que se ejecuta.
  • Lo mismo sucede para las subidas y bajadas de niveles de producción: basta con decidir para cada hora cuánto se produce, y si eso implica subir o bajar el nivel, la subida o bajada se produce al momento.
  • Justo en el cambio de una determinada hora (es decir, el el minuto 0 y segundo 0) es cuando se ejecutan las instrucciones correspondientes: encender/apagar, producir a un determinado nivel, etc. Como esas acciones se consideran instantáneas, entre el inicio de la hora $h$ y el inicio de la hora $h+1$ la situación de la central asociada no cambia.
  • Las centrales de la D a la J se considera que no tienen coste sanitario, sean cuales sean el día y la hora.

Optimización¶

A continuación vamos a ir generando el modelo de optimización correspondiente, primero con algunas simplificaciones adicionales y luego ya reproduciendo todo el funcionamiento descrito arriba.

Apartado A (4 puntos)¶

Suponiendo que todas las centrales están en marcha a las 00 horas del 1 de julio (es decir, que no hay que encenderlas en caso de que haga falta que produzcan), y que ese día tampoco se puede apagar ninguna, se pide:

  1. Formular el problema matemáticamente de cómo se debe operar a lo largo del día.
  2. Trasladar esa formulación a la librería PuLP.
  3. Resolver el problema y analizar la solución obtenida, pintando:
    • Cuánto produce cada central en cada hora.
    • ¿Hay alguna que produzca lo mínimo posible? Si es así, ¿durante qué horas del día?
    • ¿Hay alguna que produzca lo máximo posible? Si es así, ¿durante qué horas del día?
    • ¿Hay alguna que cambie entre una hora y la siguiente lo máximo posible? Si es así, ¿creciendo o decreciendo de producción?
    • Cuál es la producción total por hora.
    • ¿Cómo compara con respecto a la demanda?
    • Cuál es el coste incurrido por cada central en cada hora, desglosado en sus componentes (combustible, operación y sanitario).
    • ¿Cuál es el componente mayoritario de coste?
    • ¿Y el menos importante?
    • Cuál es el coste total por hora.
    • ¿Qué horas son las más y las menos costosas?
    • Cuál es el coste por central en el día.
    • ¿Qué centrales son las más y las menos costosas?
    • Cuál es el coste total del día.
  4. ¿Cuánto tiene que caer la demanda en % para que el problema esté en el límite de ser o no factible? Aquí puedes suponer que la demanda decae en el mismo porcentaje para las 24 horas del día. Muestra que para ese valor límite:
  • En algún momento del día la producción cae a su mínimo posible. ¿Qué momento es ese y por qué?
  • El coste es más bajo que en el escenario anterior. ¿Cuál sería este coste mínimo?
  • Si decae todavía un poco más, el problema se convierte en no factible.

1. Formulación matemática¶

1.1 Esquema visual:¶
No description has been provided for this image
Definición de variables:¶
  • $c \text{: central de producción de energía} \in \text{conjunto de centrales } C$
  • $h \text{: hora} \in \text{conjunto de horas que componen un día } H$
  • $e_{c,h}$​: Energía porducida por la central $c$, durante la hora $h$
  • $x_{c,h}$: Variable binaria que indica si la central $c$ está encendida (1) o no (0) en la hora $h$.
  • $\gamma_{c,h}$: Variable binaria que indica la existencia de un cambio de estado la central $c$ en la hora $h$ respecto a $h-1$, siendo (1) que hubo cambio de estado y (0) que no hubo.
  • $P_c$: Costo de combustible de la central $c$.
  • $O_{c}$: Costo de operativo de la central $c$ por hora.
  • $S_{c,h}$: Costo sanitario por hora de la central $c$ en la hora $h$.
  • $D_h$ Demanda total de energía en la hora $h$.
  • $ON_c$: Costo de encender/apagar de la central $c$
  • $\alpha_c$: Porcentaje de producción de la central $c$
  • $\beta_c$: Porcentaje de cambio mínimo de producción de energía por central

Donde:

Las Variables de decisión:

  • $e_{c,h}$​
  • $x_{c,h}$

Tienen la función objetivo de minimizar el costo total de producción de energía. Este costo se compone de:

  1. Costo de combustible:

    $\sum_{c\in C} \sum_{h \in H} P_c \cdot e_{c,h}$

  2. Costo operacional:

    $\sum_{c \in C} \sum_{h \in H} x_{c,h} \cdot O_{c}$

  3. Costo sanitario:

    $\sum_{c \in C} \sum_{h \in H} S_{c,h} \cdot e_{c,h}$

  4. Costo de encender o apagar una central:

    $\sum_{c \in C} \sum_{h \in H} ON_c \cdot x_{c,h}$.

La Función objetivo puede ser expresada como:

$$ \min \left( \sum_{c \in C} \sum_{h \in H} P_c \; e_{ch} + \sum_{c \in C} \sum_{h \in H} O_c \; x_{ch} + \sum_{c \in C} \sum_{h \in H} S_{ch} \; e_{ch} + \sum_{c \in C} \sum_{h \in H} ON_c \cdot \gamma_{c,h} \right) $$ $ \text{s.t.} $ $$ \sum_{c \in C} e _{c,h} \geq D_h \; \forall\; h \in H \hspace{40pt} \text{(Satisfacción de la demanda)} $$

$$ e_{c,h} \leq x_{c,h} \cdot E_c \; \forall \; h \in H, c \in H \hspace{20pt} \text{(Producción de energía máxima)} $$

$$ e_{c,h} \geq x_{c,h} \cdot E_c \cdot \alpha \; \forall \; h \in H, \hspace{30pt} \text{(Producción de energía mínima)} \\ \text{Donde: }\\ \alpha = 0.8 \; \forall \; c \in \text{Conjunto de centrales nucleares}\\ \alpha = 0.01 \; \forall \; c \notin \text{Conjunto de centrales nucleares}\\ $$

$$ \hspace{40pt} \| e_{c,h} - e_{c,h-1} \| \leq \beta \cdot E_c \; \forall \; h \in H, \hspace{20pt} \text{(Cambio máximo de producción de energía)}\\ \text{Donde: } \\ \beta = 0.2 \; \forall \; c \in \text{Conjunto de centrales nucleares} \\ \beta = 0.25 \; \forall \; c \in \text{Conjunto de centrales a carbón} \\ \beta = 1 \; \forall \; c \notin \text{Conjunto de centrales nucleares o carbón} $$

$$ x_{c,h} = 1 \; \forall \; h \in H , c \in C\;\hspace{10pt} \text{(Estado de plantas)}\\ $$

$$ \gamma_{c,h} = \| x_{c,h-1} - x_{c,h} \| \hspace{15pt} \text{(Cambio de estado de encendido/apagado)}\\ $$

Ahora bien, la función objetivo se define de manera general para todo el ejercicio. Sin embargo, debido a que se supone que todas las centrales están en marcha a las 00 horas del 1 de julio (es decir, que no hay que encenderlas en caso de que haga falta que produzcan), y que ese día tampoco se puede apagar ninguna, la variable $x_{c,h}$ siempre es 1.

2. Formulación en la librería PulP¶

In [11]:
# Definición del problema
prob = LpProblem("Optimizacion_Energia", LpMinimize)
Definición de variables¶
In [12]:
C = list(plants.reset_index()['central']) # Conjunto de centrales C
H = list(daily_demand.index) # Conjunto de horas que componen un día H
In [13]:
def crear_variables_decision(C, H):
    # Variables de decisión
    e_ch = LpVariable.dicts('MWh', (C, H), lowBound=0, cat=LpContinuous) # Energía producida por la central c durante la hora h
    x_ch = LpVariable.dicts('X', (C, H), cat='Binary') # Variable binaria que indica si la central c esta encendida o no en la hora h
    gamma_ch = LpVariable.dicts('gamma', (C, H), cat='Binary') # Variable binaria que indica si la central c está encendida (1) o no (0) en la hora h.
    return e_ch, x_ch, gamma_ch

e_ch, x_ch, gamma_ch = crear_variables_decision(C, H)
In [14]:
alpha_values = [0.8 if c in ['D', 'E'] else 0.01 for c in C]
df_alpha = pd.DataFrame({
    'central': C,
    'alpha': alpha_values
}).set_index('central')

alpha = df_alpha.squeeze() #Porcentaje de producción de la central c


beta_values = [0.2 if c in ['D', 'E'] else 0.25 if c in ['A','B', 'C', 'J'] else 1 for c in C]
df_beta = pd.DataFrame({
    'central': C,
    'beta': beta_values
    }).set_index('central')

beta = df_beta.squeeze() #Porcentaje de cambio mínimo de producción de energía por central
In [15]:
#Producción máxima
E_c = plants.reset_index()[['central','producción']].set_index('central') #Energía máxima de energía para una central c
In [16]:
## Costo del combustible
P_c = plants.join(plant_costs, on= 'combustible').reset_index()[['central','coste_comb']].set_index('central') # Unimos las tablas y seleccionamos las columnas necesarias
P_c.index.name = None # Unimos las tablas y seleccionamos las columnas necesarias
for h in H:
    P_c[h] = P_c['coste_comb']
P_c = P_c.drop(['coste_comb'], axis=1) # Costo de combustible de la central $c$.
In [17]:
## Costos operativos
O_c = plants.join(plant_costs, on= 'combustible').reset_index()[['central','coste_oper']].set_index('central')
O_c.index.name = None # Unimos las tablas y seleccionamos las columnas necesarias
for h in H:
    O_c[h] = O_c['coste_oper']
O_c = O_c.drop(['coste_oper'], axis=1) # Costo de operativo de la central $c$ por hora.
In [18]:
## Costos sanitarios 
S_ch = health_costs.reset_index().sort_index() # Costo sanitario por hora de la central $c$ en la hora $h$.
S_ch = S_ch.pivot(index='hora', columns='central', values ='coste_salud') # Pivotamos el DataFrame 'S_ch' para reorganizarlo y facilitar su análisis.

non_carbon_plants = ['D', 'E', 'F', 'G', 'H', 'I', 'J']

# El método reindex se utiliza para agregar las plantas sin costos de salud al DataFrame.
# El parámetro 'fill_value' se establece en 0, lo que reemplaza todos los valores faltantes con 0.
S_ch= S_ch.reindex(columns=S_ch.columns.tolist() + non_carbon_plants, fill_value=0)
In [19]:
D_h = daily_demand # Demanda total de energía en la hora $h$.
In [20]:
ON_c = plants.join(plant_costs, on= 'combustible').reset_index()[['central','coste_on_off']].set_index('central')
ON_c.index.name = None # Unimos las tablas y seleccionamos las columnas necesarias
for n in H:
    ON_c[n] = ON_c['coste_on_off']
ON_c = ON_c.drop(['coste_on_off'], axis=1) # Costo de encender/apagar de la central $c$
Definición de la función objetivo¶
In [21]:
def crear_funcion_obj(prob, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, C, H):
    # Función: Costo de combustible
    comb = (lpSum([P_c[h][c] * e_ch[c][h] for c in C for h in H]))

#Función : Costo operacional
    op =  (lpSum([O_c[h][c] * x_ch[c][h] for c in C for h in H]))

# Función: Costo sanitario
    san = lpSum([S_ch[c][h] * e_ch[c][h] for c in C for h in H])

# Función costo de encender o apagar una central (Se define para su uso a posteriori, puesto que para este ejercicio nunca se aplicará)

    on_off = (lpSum([ON_c[h][c] * gamma_ch[c][h] for c in C for h in H]))

# Definiendo la función objetivo como:
    
    prob+= comb + op + san + on_off
    
    return prob 
In [22]:
prob= crear_funcion_obj(prob, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, C, H)
Definición de las restricciones¶
In [23]:
def anadir_restriccion_demanda(prob, e_ch, D_h, C, H):
    ## Satisfacción de la demanda
    for h in H:
        prob += (lpSum(e_ch[c][h] for c in C) >= D_h.loc[h][0]) # Satisfacción de la demanda

    return prob
In [24]:
prob= anadir_restriccion_demanda(prob, e_ch, D_h, C, H)
In [25]:
def anadir_restriccion_xch(prob, x_ch, C, H):
    # Todas las centrales están encendidas en todas las horas
    for c in C:
        for h in H:
            prob += x_ch[c][h] == 1  
    
    return prob
In [26]:
prob= anadir_restriccion_xch(prob, x_ch, C, H)
In [27]:
def anadir_otras_restricciones(prob, H, C, e_ch, x_ch, E_c, beta, alpha, gamma_ch):
    ## Producción máxima
    for c in C:
        for h in H:
            prob += e_ch[c][h] <= x_ch[c][h] * E_c.loc[c][0]

    ## Producción mínima
    for c in C:
        for h in H:
            prob += e_ch[c][h] >= x_ch[c][h] * E_c.loc[c][0]*alpha[c]

    ## Cambio máximo de producción de energia
    for c in C:
        for h in range(1, len(H[0:-1])):
            prob += e_ch[c][H[h]] - e_ch[c][H[h-1]] <= beta[c] * E_c.loc[c][0]
            prob += e_ch[c][H[h-1]] - e_ch[c][H[h]] <= beta[c] * E_c.loc[c][0]    

    # Cambio de estado de cada central
    for c in C:
        for h in range(1, len(H[0:-1])):
            prob += gamma_ch[c][H[h]] == x_ch[c][H[h]] - x_ch[c][H[h-1]]  # Cambio de estado de encendido/apagado
            prob += gamma_ch[c][H[h]] == x_ch[c][H[h-1]] - x_ch[c][H[h]]  # Cambio de estado de encendido/apagado           

    return prob                 
In [28]:
prob= anadir_otras_restricciones(prob, H, C, e_ch, x_ch, E_c, beta, alpha, gamma_ch)
In [29]:
status = prob.solve()
print('Status:', LpStatus[status])
value(prob.objective)
costo_inicial = value(prob.objective)
Status: Optimal

3. Resolver el problema y analizar la solución obtenida, pintando:¶

In [30]:
def almacenar_resultados_modelo(prob):   
    energy = {}
    x = {}
    gamma = {}
    result_plant = []
    result_energy = []
    result_hour = []
    block_size = 240
    block_count = 0
    current_dict = energy  # Start with energy

    for v in prob.variables():
        p = (v.name).split('_')[1]
        e = v.varValue
        h = (v.name).split('_')[2] + '-' + (v.name).split('_')[3] + '-' + (v.name).split('_')[4] + '-' + (v.name).split('_')[5]

        result_plant.append(p)
        result_energy.append(e)
        result_hour.append(h)

        block_count += 1
        if block_count % block_size == 0:
            if current_dict == energy:
                current_dict = x
            elif current_dict == x:
                current_dict = gamma
            else:
                break

            result_plant = []  
            result_energy = []
            result_hour = []

        current_dict['central'] = result_plant
        current_dict['hora'] = result_hour
        current_dict['energía'] = result_energy

    df_energy = pd.DataFrame(energy)
    df_xch = pd.DataFrame(x)
    df_gamma = pd.DataFrame(gamma)

    return (df_energy, df_xch, df_gamma)
In [31]:
df_result_energy, df_result_xch, df__result_gamma= almacenar_resultados_modelo(prob)

Cuánto produce cada central en cada hora.¶
In [32]:
centrals = df_result_energy['central'].unique()

plt.figure(figsize=(20,10))

for central in centrals:
    plt.plot(df_result_energy[df_result_energy['central'] == central]['hora'], 
             df_result_energy[df_result_energy['central'] == central]['energía'], 
             marker='o', linestyle='-', 
             label=f'Energía producida por la central {central}')

    
plt.title(f'Energía producida por cada central por hora', weight= 'bold')
plt.xlabel('Hora')
plt.ylabel('Energía producida (MWh)')
plt.tick_params(axis='x', rotation=90)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()
No description has been provided for this image
¿Hay alguna que produzca lo mínimo posible? Si es así, ¿durante qué horas del día?, ¿Hay alguna que produzca lo máximo posible? Si es así, ¿durante qué horas del día?¶
In [33]:
centrals = df_result_energy['central'].unique()
fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2
    limite_superior = E_c.loc[central, 'producción']  # Obtiene el límite superior para la central
    limite_inferior = E_c.loc[central, 'producción'] * df_alpha.loc[central, 'alpha']  # Obtiene el límite inferior para la central
    axs[row, col].plot(df_result_energy[df_result_energy['central'] == central]['hora'], 
                       df_result_energy[df_result_energy['central'] == central]['energía'], 
                       marker='o', linestyle='-', 
                       label=f'Energía producida por la central {central}')
    axs[row, col].axhline(y=limite_superior, color='r', linestyle='--', label='Capacidad máxima de la central')
    axs[row, col].axhline(y=limite_inferior, color='b', linestyle='--', label='Mínima produccción admitida')
    axs[row, col].set_title(f'Producción de energia por central {central}')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Energía producida (MWh)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()
No description has been provided for this image
¿Hay alguna que cambie entre una hora y la siguiente lo máximo posible? Si es así, ¿creciendo o decreciendo de producción?¶
In [34]:
fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2
    df_central = df_result_energy[df_result_energy['central'] == central]
    df_central['cambio_energia'] = df_central['energía'].diff()  # Calcula la derivada de la energía
    
    limite_superior = E_c.loc[central, 'producción'] * df_beta.loc[central, 'beta']  # Obtiene el límite superior para la central
    limite_inferior = E_c.loc[central, 'producción'] * df_beta.loc[central, 'beta']  * -1 # Obtiene el límite inferior para la central
    
    axs[row, col].axhline(y=limite_superior, color='r', linestyle='--', label='Incremento de producción máximo admisible de la central')
    axs[row, col].axhline(y=limite_inferior, color='b', linestyle='--', label='Descenso de producción máximo admisible de la central')

    axs[row, col].plot(df_central['hora'], 
                       df_central['cambio_energia'], 
                       marker='o', linestyle='-', 
                       label=f'Cambio de energía en la central {central}')
    axs[row, col].set_title(f'Cambio de producción de energía por hora en la central {central}')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Cambio de energía producida (MWh/h)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()
No description has been provided for this image

Por medio de la derivada de la funcion de producción de energía es posible visualizar los incrementos y descensos realizados por cada planta en su produccicón; valores en cero indican que la producción fue constante entre una hora y la siguiente. A su vez se establecen para cada uno de los gráficos los limites superiores e inferiores admisibles correspondientes a cada planta de acuerdo a su tipología de combustible: Los graficos demuestran que en ninguno de los casos son sobrepasados o igualados (pasar de 0 a alguno de los limites) los limites de incremento y descenso de producción de energía admisibles.

Cuál es la producción total por hora¶
In [35]:
df_result_energy['hora'] = pd.to_datetime(df_result_energy['hora'])
df_energia_demanda = pd.merge(df_result_energy, D_h, on='hora', how='left')
df_energia_demanda = df_energia_demanda.astype({'hora': 'string'})
In [36]:
produccion_total = df_energia_demanda.groupby(df_energia_demanda['hora'])[['energía', 'demanda']].sum().reset_index()
produccion_total['demanda'] = produccion_total['demanda']/10

plt.figure(figsize=(18,8))

plt.plot(produccion_total['hora'], produccion_total['energía'], label = 'Producción de energía total por hora', marker='x')

plt.xlabel('Hora')
plt.ylabel('MWh')
plt.xticks(rotation=90)

plt.legend()
plt.show
Out[36]:
<function matplotlib.pyplot.show(close=None, block=None)>
No description has been provided for this image
¿Cómo compara con respecto a la demanda?¶
In [37]:
plt.figure(figsize=(18,8))

plt.plot(produccion_total['hora'], produccion_total['energía'], label = 'Producción de energía total por hora', marker='x')
plt.plot(produccion_total['hora'], produccion_total['demanda'], label = 'Demanda de energía total por hora', linestyle = '--')

plt.fill_between(produccion_total['hora'], produccion_total['energía'], produccion_total['demanda'], 
                 where = (produccion_total['energía'] >= produccion_total['demanda']), facecolor='green', alpha = 0.5)

plt.xlabel('Hora')
plt.ylabel('MWh')
plt.xticks(rotation=90)

plt.legend()
plt.show
Out[37]:
<function matplotlib.pyplot.show(close=None, block=None)>
No description has been provided for this image

Hay una coincidencia casi perfecta. Pero al mirar con lupa los datos en formato tabla, se evidencian diferencias. Visualmente estas diferencias son:

In [38]:
# Calcular la diferencia entre la producción y la demanda
produccion_total['diferencia'] = produccion_total['energía'] - produccion_total['demanda']

plt.figure(figsize=(10,5))

plt.plot(produccion_total['hora'], produccion_total['diferencia'], label='Diferencia', color='purple')

plt.title('Diferencia entre la producción y la demanda')
plt.xlabel('Hora')
plt.ylabel('Diferencia (MWh)')
plt.legend()
plt.xticks(rotation=90)


plt.show()
No description has been provided for this image
Cuál es el coste incurrido por cada central en cada hora, desglosado en sus componentes (combustible, operación y sanitario).¶
In [39]:
# Transformamos el dataframe P_c a un formato largo
P_c_long = P_c.reset_index().melt(id_vars='index', var_name='hora', value_name='costo_combustible')
P_c_long.rename(columns={'index':'central'}, inplace=True)
P_c_long = P_c_long.astype({'hora': 'string'})
In [40]:
costos = pd.merge(df_energia_demanda, P_c_long,  how='left', left_on=['central','hora'], right_on = ['central','hora'])
costos['costo_combustible'] = costos['energía'] * costos ['costo_combustible']
In [41]:
O_c_long = O_c.reset_index().melt(id_vars='index', var_name='hora', value_name='costo_operativo')
O_c_long.rename(columns={'index':'central'}, inplace=True)
O_c_long = O_c_long.astype({'hora': 'string'})
In [42]:
costos = pd.merge(costos, O_c_long,  how='left', left_on=['central','hora'], right_on = ['central','hora'])
In [43]:
S_ch_T = S_ch.T
S_ch_T_long = S_ch_T.reset_index().melt(id_vars='central', var_name='hora', value_name='costo_sanitario')
S_ch_T_long.rename(columns={'index':'central'}, inplace=True)
S_ch_T_long = S_ch_T_long.astype({'hora': 'string'})
In [44]:
costos = pd.merge(costos, S_ch_T_long,  how='left', left_on=['central','hora'], right_on = ['central','hora'])
costos['costo_sanitario'] = costos['energía'] * costos ['costo_sanitario']
In [45]:
fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2
    df_central = costos[costos['central'] == central]

    df_central.groupby('hora').sum()[['costo_operativo', 'costo_combustible', 'costo_sanitario']].plot(kind='bar', stacked=True, ax=axs[row, col])
    
    axs[row, col].set_title(f'Distribución de costos de produccion de energía por hora en la central {central}')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Costo (Costo/MWh)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()
No description has been provided for this image
Cuál es el coste total por hora.¶
In [46]:
costos['costo total'] = costos['costo_operativo'] + costos['costo_combustible'] + costos['costo_sanitario']
costo_total = costos.groupby(costos['hora'])['costo total'].sum().reset_index()
In [47]:
# Crea un gráfico de barras con un ancho adecuado
plt.figure(figsize=(20, 6))
plt.bar(costo_total['hora'], costo_total['costo total'], label ='Costo total por hora')

plt.xlabel('Horas')
plt.ylabel('Costos')
plt.title('Costo Total por hora')
plt.legend()
plt.xticks(rotation=90)
plt.show()
No description has been provided for this image
¿Qué horas son las más y las menos costosas?¶
In [48]:
costo_minimo = costo_total['costo total'].min()
costo_maximo = costo_total['costo total'].max()
hora_minima = costo_total['hora'].loc[costo_total['costo total']==costo_minimo].values[0]
hora_maxima = costo_total['hora'].loc[costo_total['costo total']==costo_maximo].values[0]
In [49]:
print(f'La hora menos costosa es: {hora_minima}')
La hora menos costosa es: 2023-07-01 09:00:00
In [50]:
print(f'La hora más costosa es: {hora_maxima}')
La hora más costosa es: 2023-07-01 20:00:00
Cuál es el coste por central en el día.¶
In [51]:
costo_total_por_planta = costos.groupby(costos['central'])['costo total'].sum().reset_index()
plt.figure(figsize=(20, 6))
plt.bar(costo_total_por_planta['central'], costo_total_por_planta['costo total'], label ='Costo total por central')

plt.xlabel('Centrales')
plt.ylabel('Costos')
plt.title('Costo Total por central al día')
plt.legend()
plt.show()
No description has been provided for this image
In [52]:
costo_minimo = costo_total_por_planta['costo total'].min()
costo_maximo = costo_total_por_planta['costo total'].max()
central_min = costo_total_por_planta['central'].loc[costo_total_por_planta['costo total']==costo_minimo].values[0]
central_max = costo_total_por_planta['central'].loc[costo_total_por_planta['costo total']==costo_maximo].values[0]
¿Qué centrales son las más y las menos costosas?¶
In [53]:
print(f'La central menos costosa es: {central_min}')
La central menos costosa es: I
In [54]:
print(f'La central más costosa es: {central_max}')
La central más costosa es: A
Cuál es el coste total del día¶
In [55]:
coste_dia = value(prob.objective)
print(f'El coste total del día: {coste_dia}')
El coste total del día: 4508697.3013742585

4. ¿Cuánto tiene que caer la demanda en % para que el problema esté en el límite de ser o no factible? Aquí puedes suponer que la demanda decae en el mismo porcentaje para las 24 horas del día. Muestra que para ese valor límite:¶

In [56]:
def hallar_limite_viable(demanda, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha):
    prob = LpProblem('Problema_de_la_demanda', LpMinimize)
    e_ch, x_ch, gamma_ch = crear_variables_decision(C, H)
    prob = crear_funcion_obj(prob, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, C, H)
    prob = anadir_restriccion_xch(prob, x_ch, C, H)
    prob = anadir_otras_restricciones(prob, H, C, e_ch, x_ch, E_c, beta, alpha, gamma_ch)
    for h in H:
        prob+= (lpSum([e_ch[c][h] for c in C])) == demanda.loc[h][0]

    return prob, e_ch, x_ch, gamma_ch
In [57]:
coeficiente = 1
paso = 0.01
status = 'Optimal'

while (status == 'Optimal'):
    demanda = daily_demand * coeficiente
    prob, e_ch, x_ch, gamma_ch = hallar_limite_viable(demanda, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha)
    status = LpStatus[prob.solve()]
    coeficiente-= paso    
In [58]:
porcentaje_inviable = 100 - (coeficiente+paso)*100
print(f"El porcentaje de reducción de demanda de energía que hace que la solución sea inviable es: ", porcentaje_inviable, "%")
El porcentaje de reducción de demanda de energía que hace que la solución sea inviable es:  58.00000000000005 %
In [59]:
limite_viable = daily_demand * (coeficiente + 2*paso)
prob, e_ch, x_ch, gamma_ch = hallar_limite_viable(limite_viable, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha)
In [60]:
status = prob.solve()
print('Status:', LpStatus[status])
value(prob.objective)
Status: Optimal
Out[60]:
910405.7955026062
In [61]:
df_result_limits, df_result_xch, df_result_gamma = almacenar_resultados_modelo(prob)
df_result_limits['hora'] = pd.to_datetime(df_result_limits['hora'])
df_energia_demanda_limite = pd.merge(df_result_limits, D_h, on='hora', how='left')
df_energia_demanda_limite = df_energia_demanda_limite.astype({'hora': 'string'})
En algún momento del día la producción cae a su mínimo posible. ¿Qué momento es ese y por qué?¶
In [62]:
produccion_total_limite = df_energia_demanda_limite.groupby(df_energia_demanda_limite['hora'])[['energía', 'demanda']].sum().reset_index()
produccion_total_limite['demanda'] = produccion_total_limite['demanda']/10

plt.figure(figsize=(18,8))

plt.plot(produccion_total_limite['hora'], produccion_total_limite['energía'], label = 'Producción de energía total por hora', marker='x')

plt.xlabel('Hora')
plt.ylabel('MWh')
plt.xticks(rotation=90)

plt.legend()
plt.show
Out[62]:
<function matplotlib.pyplot.show(close=None, block=None)>
No description has been provided for this image
In [63]:
produccion_minima = produccion_total_limite['energía'].min()
hora_minima = produccion_total_limite['hora'].loc[produccion_total_limite['energía']==produccion_minima].values[0]
print(f'El momento en que la producción cae a su mínimo posible es: {hora_minima}')
El momento en que la producción cae a su mínimo posible es: 2023-07-01 08:00:00
In [64]:
demanda_minima = produccion_total_limite['demanda'].min()
hora_minima = produccion_total_limite['hora'].loc[produccion_total_limite['demanda']==demanda_minima].values[0]
print(f'El momento en que la demanda es menor es: {hora_minima}')
El momento en que la demanda es menor es: 2023-07-01 08:00:00

El momento del día en que la producción cae a su mínimo posible es a las 8 AM, debido a que este es el momento con menor demanda energética.

El coste es más bajo que en el escenario anterior. ¿Cuál sería este coste mínimo?¶
In [65]:
status = prob.solve()
print("Status: ", LpStatus[status])
costo_limite = value(prob.objective) 
value(prob.objective)
Status:  Optimal
Out[65]:
910405.7955026062
In [66]:
dif =round((100 - ((costo_limite*100)/costo_inicial)),2)
print(dif)
79.81
In [67]:
print(f'El costo mínimo para este escenario es {costo_limite}, comparado con el costo inicial de {costo_inicial}, lo que indica una caída del {dif}% en el costo')
El costo mínimo para este escenario es 910405.7955026062, comparado con el costo inicial de 4508697.3013742585, lo que indica una caída del 79.81% en el costo
Si decae todavía un poco más, el problema se convierte en no factible.¶
In [68]:
limite_viable = daily_demand * (coeficiente + paso)
prob, e_ch, x_ch, gamma_ch = hallar_limite_viable(limite_viable, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha)
In [69]:
status = prob.solve()
print('Status:', LpStatus[status])
value(prob.objective)
Status: Infeasible
Out[69]:
875612.3748042658

Apartado B (6 puntos)¶

Volviendo a la demanda original, y sin suponer ya que las centrales no se puedan apagar (excepto las nucleares, que por restricción tienen que operar siempre), ni tampoco que al principio tengan por qué estar todas encendidas, se pide:

  1. Formular el problema matemáticamente de cómo se debe operar a lo largo del día.
  • ¿Qué cambios hace falta introducir con respecto a la formulación del apartado A?
  1. Trasladar esta nueva formulación a la librería PuLP.
  2. Resolver el problema y analizar la solución obtenida, comparando con el resultado obtenido antes:
    • ¿Cómo varía la producción de energía de cada central? La demanda sigue siendo la misma, pero es posible que no se opere exactamente igual que antes...
    • ¿Hay alguna central que se desactive o active? Si es así, ¿en qué horas?
    • ¿Hay alguna siempre encendida (además de las nucleares)?
    • ¿Y alguna siempre apagada?
    • ¿Cómo varía la distribución de costes? ¿Es más barato o más caro que antes el total del día? ¿Por qué?
  3. ¿Cuánto tiene que subir la demanda en % para que el problema esté en el límite de ser o no factible? De nuevo supón que la demanda aumenta en el mismo porcentaje para las 24 horas del día. Muestra que para ese valor límite:
  • En algún momento del día la producción sube al máximo posible. ¿Qué momento es ese y por qué?
  • El coste es mayor que con la demanda predicha. ¿Cuál sería ese coste máximo?
  • Si sube todavía un poco más, el problema se convierte en no factible.

1. Formular el problema matemáticamente de cómo se debe operar a lo largo del dia¶

La Función objetivo puede ser expresada como:

$$ \min \left( \sum_{c \in C} \sum_{h \in H} P_c \; e_{ch} + \sum_{c \in C} \sum_{h \in H} O_c \; x_{ch} + \sum_{c \in C} \sum_{h \in H} S_{ch} \; e_{ch} + \sum_{c \in C} \sum_{h \in H} ON_c \cdot \gamma_{c,h} \right) $$ $ \text{s.t.} $ $$ \sum_{c \in C} e _{c,h} \geq D_h \; \forall\; h \in H \hspace{35pt} \text{(Satisfacción de la demanda)} $$

$$ e_{c,h} \leq x_{c,h} \cdot E_c \; \forall \; h \in H, c \in H \hspace{15pt} \text{(Producción de energía máxima)} $$

$$ e_{c,h} \leq x_{c,h} \cdot E_c \cdot \alpha \; \forall \; h \in H, \hspace{25pt} \text{(Producción de energía mínima)} \\ \text{Donde: }\\ \alpha = 0.8 \; \forall \; c \in \text{Conjunto de centrales nucleares}\\ \alpha = 0.01 \; \forall \; c \notin \text{Conjunto de centrales nucleares}\\ $$

$$ \hspace{50pt} \| e_{c,h} - e_{c,h-1} \| \leq \beta \cdot E_c \; \forall \; h \in H, \hspace{15pt} \text{(Cambio máximo de producción de energía)}\\ \text{Donde: } \\ \beta = 0.2 \; \forall \; c \in \text{Conjunto de centrales nucleares} \\ \beta = 0.25 \; \forall \; c \in \text{Conjunto de centrales a carbón} \\ \beta = 1 \; \forall \; c \notin \text{Conjunto de centrales nucleares o carbón} $$

$$ x_{c,h} = 1 \; \forall \; h \in H , c \in \text{Conjunto de centrales nucleares} \;\hspace{15pt} \text{(Estado de plantas nucleares)}\\ $$

$$ \gamma_{c,h} = \| x_{c,h-1} - x_{c,h} \| \hspace{15pt} \text{(Cambio de estado de encendido/apagado)}\\ $$

  • ¿Qué cambios hace falta introducir con respecto a la formulación del apartado A?

Debido a que en el apartado 1 se definió la función de manera general, no se deben realizar cambios al planteamiento de la formulación matemática. En este caso particular, a la hora de realizar la implementación, no se definirá Xc,h como 1 para todas las centrales, debido a que pueden existir cambios de estados, ni se asume que todas inician encendidas.

2. Trasladar esta nueva formulación a la librería PuLP.¶

In [70]:
# Definicion del problema
prob = LpProblem('Optimizacion_Energia_ONOFF', LpMinimize)
Variables de decisión¶
In [71]:
# Debido a que el codigo del apartado anterior se definió modularmente y general, entonces:
#Creamos las variables decisión
e_ch, x_ch, gamma_ch = crear_variables_decision(C, H)

#No es necesario redifinir las variables, puesto que se mantiene el contexto.
Función Objetivo¶
In [72]:
# Definicion de la función objetivo:
prob= crear_funcion_obj(prob, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, C, H)
Definición de las restricciones¶

Este es el único apartado que es diferente a la formulación anterior, puesto que la restricción relacionada a la variable x_ch cambia. Por tanto, no se llama a su función y se define una nueva.

In [73]:
prob= anadir_restriccion_demanda(prob, e_ch, D_h, C, H) # Satisfacción de la demanda
In [74]:
def anadir_restriccion_xch2(prob, x_ch, C, H):
    centrales_nucleares = ['D', 'E']
    for c in C:
        for h in H:
            if (c in centrales_nucleares):
                prob += x_ch[c][h]==1 # Las centrales nucleares siempre deben estar encendidas
    return prob
In [75]:
prob = anadir_restriccion_xch2(prob, x_ch, C, H) 
In [76]:
prob= anadir_otras_restricciones(prob, H, C, e_ch, x_ch, E_c, beta, alpha, gamma_ch)
In [77]:
status = prob.solve()
print('Status:', LpStatus[status])
value(prob.objective)
costo_inicial_onoff = value(prob.objective)
costo_inicial_onoff
Status: Optimal
Out[77]:
4497165.960429066

3. Resolver el problema y analizar la solución obtenida, comparando con el resultado obtenido antes:¶

In [78]:
df_results, df_result_xch, df_result_gamma = almacenar_resultados_modelo(prob)
¿Cómo varía la producción de energía de cada central? La demanda sigue siendo la misma, pero es posible que no se opere exactamente igual que antes..¶
In [79]:
centrals = df_results['central'].unique()

plt.figure(figsize=(20,10))

for central in centrals:
    plt.plot(df_results[df_results['central'] == central]['hora'], 
             df_results[df_results['central'] == central]['energía'], 
             marker='o', linestyle='-', 
             label=f'Energía producida por la central {central}')

    
plt.title(f'Energía producida por cada central por hora', weight= 'bold')
plt.xlabel('Hora')
plt.ylabel('Energía producida (MWh)')
plt.tick_params(axis='x', rotation=90)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.tight_layout()
plt.show()
No description has been provided for this image

A simple vista, se pueden apreciar cambios en la producción de la central B. Sin embargo, vale la pena revisarlo con lupa. Por tanto:

In [80]:
df_results = df_results.astype({'hora': 'datetime64[ns]'})
compared_results = pd.merge(df_results, df_result_energy, how= 'left', on=['central', 'hora'])
In [81]:
centrals = df_results['central'].unique()
fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2
    limite_superior = E_c.loc[central, 'producción']  # Obtiene el límite superior para la central
    limite_inferior = E_c.loc[central, 'producción'] * df_alpha.loc[central, 'alpha']  # Obtiene el límite inferior para la central
    axs[row, col].plot(df_results[df_results['central'] == central]['hora'], 
                       df_results[df_results['central'] == central]['energía'], 
                       marker='o', linestyle='-', 
                       label=f'Energía producida por la central {central}')
    axs[row, col].axhline(y=limite_superior, color='r', linestyle='--', label='Capacidad máxima de la central')
    axs[row, col].axhline(y=limite_inferior, color='b', linestyle='--', label='Mínima produccción admitida')
    axs[row, col].set_title(f'Producción de energia por central {central}')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Energía producida (MWh)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()
No description has been provided for this image
In [82]:
fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2
    df_central = compared_results[compared_results['central'] == central]

    axs[row, col].plot(df_central['hora'], 
                       df_central['energía_x'], 
                       marker='o', linestyle='-', 
                       label=f'Producción de energia con centrales encendidas y apagadas')
    
    axs[row, col].plot(df_central['hora'], 
                       df_central['energía_y'], 
                       marker='o', linestyle='--', 
                       label=f'Producción de energia con centrales siempre encendidas')
    
    axs[row, col].fill_between(df_central['hora'], df_central['energía_y'], df_central['energía_x'], 
                               where = (df_central['energía_y'] != df_central['energía_x']), facecolor= 'gray', alpha=0.5 )

    axs[row, col].set_title(f'Diferencia de energía producida en la central {central}, respecto al escenario anterior')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Energía producida (MWh/h)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()
No description has been provided for this image

Se puede observar, entonces, que hay variaciones muy pequeñas entre las producciones de las centrales A, B, C, H y J. Las centrales D, E, F, G e I mantienen la misma producción relativa al ejercicio anterior.

¿Hay alguna central que se desactive o active? Si es así, ¿en qué horas?¶
In [83]:
centrals = df_result_xch['central'].unique()
fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2

    axs[row, col].plot(df_result_xch[df_result_xch['central'] == central]['hora'], 
                       df_result_xch[df_result_xch['central'] == central]['energía'], 
                       marker='o', linestyle='-', 
                       label=f'Energía producida por la central {central}')
    axs[row, col].axhline(y=1, color='r', linestyle='--', label='Capacidad máxima de la central')
    axs[row, col].axhline(y=0, color='b', linestyle='--', label='Mínima produccción admitida')
    axs[row, col].set_title(f'Producción de energia por central {central}')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Energía producida (MWh)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()
No description has been provided for this image

Los resultados de la solución óptima refleja que las centrales C y B se apagan desde un principio, y no vuelven a ser encendidas. Para el resto de centrales, se evidencia que inician encendidas y permanecen en este estado durante todo el día

¿Hay alguna siempre encendida (además de las nucleares)?¶

Además de las nucleares, se cuentran siempre encendidas las centrales A, F, G, H, I y J

¿Y alguna siempre apagada?¶

Las centrales que permanecen siempre apagadas son la B y la C

¿Cómo varía la distribución de costes? ¿Es más barato o más caro que antes el total del día? ¿Por qué?¶
In [84]:
df_results = df_results.astype({'hora': 'string'})
In [85]:
new_costos = pd.merge(df_results, P_c_long, how='left', left_on=['central','hora'], right_on = ['central','hora'])
new_costos['costo_combustible'] = new_costos['energía'] * new_costos ['costo_combustible']
In [86]:
df_result_xch = df_result_xch.astype({'hora': 'datetime64[ns]'})
df_result_xch = df_result_xch.astype({'hora': 'string'})
In [87]:
new_costos = pd.merge(new_costos, O_c_long,  how='left', left_on=['central','hora'], right_on = ['central','hora'])
new_costos = pd.merge(new_costos, df_result_xch, how='left', on=['central', 'hora'])
In [88]:
new_costos['costo_operativo'] = new_costos['costo_operativo'] * new_costos['energía_y']
In [89]:
new_costos = pd.merge(new_costos, S_ch_T_long,  how='left', left_on=['central','hora'], right_on = ['central','hora'])
new_costos['costo_sanitario'] = new_costos['energía_x'] * new_costos ['costo_sanitario']
In [90]:
fig, axs = plt.subplots(nrows=5, ncols=2, figsize=(20,30))

for i, central in enumerate(centrals):
    row = i // 2
    col = i % 2
    df_central = new_costos[new_costos['central'] == central]

    df_central.groupby('hora').sum()[['costo_operativo', 'costo_combustible', 'costo_sanitario']].plot(kind='bar', stacked=True, ax=axs[row, col])
    
    axs[row, col].set_title(f'Distribución de costos de producción de energía por hora para la central {central}')
    axs[row, col].set_xlabel('Hora')
    axs[row, col].set_ylabel('Costo (Costo/MWh)')
    axs[row, col].tick_params(axis='x', rotation=90)
    axs[row, col].legend()

plt.tight_layout()
plt.show()
No description has been provided for this image
In [91]:
new_costos['costo_total'] = new_costos['costo_operativo'] + new_costos['costo_combustible'] + new_costos['costo_sanitario']
new_costo_total = new_costos.groupby(new_costos['hora'])['costo_total'].sum().reset_index()
In [92]:
plt.figure(figsize=(20, 6))
plt.plot(costo_total['hora'], costo_total['costo total'], label ='Costo total por hora, anterior')
plt.plot(new_costo_total['hora'], new_costo_total['costo_total'], label ='Costo total por hora actual')

plt.xlabel('Horas')
plt.ylabel('Costos')
plt.title('Costo Total por hora')
plt.legend()
plt.xticks(rotation=90)
plt.show()
No description has been provided for this image
In [93]:
new_costo_total['diferencia'] = costo_total['costo total'] - new_costo_total['costo_total']

plt.figure(figsize=(10,5))

plt.plot(new_costo_total['hora'], new_costo_total['diferencia'], label='Diferencia', color='purple')

plt.title('Diferencia entre la producción y la demanda')
plt.xlabel('Hora')
plt.ylabel('Diferencia (MWh)')
plt.legend()
plt.xticks(rotation=90)


plt.show()
No description has been provided for this image

Los reusltados reflejan que el nuevo modelo tiene una reducción de costos totales resapecto al modelo anterior, especialmente en el rango temporal entre las 15:00 y las 00:00

4. ¿Cuánto tiene que subir la demanda en % para que el problema esté en el límite de ser o no factible? De nuevo supón que la demanda aumenta en el mismo porcentaje para las 24 horas del día. Muestra que para ese valor límite:¶

In [94]:
def hallar_limite_viable(demanda, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha):
    prob = LpProblem('Problema_de_la_demanda', LpMinimize)
    e_ch, x_ch, gamma_ch = crear_variables_decision(C, H)
    prob = crear_funcion_obj(prob, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, C, H)
    prob = anadir_restriccion_xch2(prob, x_ch, C, H)
    prob = anadir_otras_restricciones(prob, H, C, e_ch, x_ch, E_c, beta, alpha, gamma_ch)
    for h in H:
        prob+= (lpSum([e_ch[c][h] for c in C])) == demanda.loc[h][0]

    return prob, e_ch, x_ch, gamma_ch
In [95]:
coeficiente = 1
paso = 0.005
status = 'Optimal'

while (status == 'Optimal'):
    demanda = daily_demand * coeficiente
    prob, e_ch, x_ch, gamma_ch = hallar_limite_viable(demanda, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha)
    status = LpStatus[prob.solve()]
    coeficiente+= paso    
In [96]:
porcentaje_inviable = ((coeficiente+paso)*100)-100
print(f"El porcentaje de reducción de demanda de energía que hace que la solución sea inviable es: ", porcentaje_inviable, "%")
El porcentaje de reducción de demanda de energía que hace que la solución sea inviable es:  32.99999999999929 %
In [97]:
limite_viable = daily_demand * 1.315
prob, e_ch, x_ch, gamma_ch = hallar_limite_viable(limite_viable, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha)
In [98]:
status = prob.solve()
print('Status:', LpStatus[status])
value(prob.objective)
Status: Optimal
Out[98]:
7755835.807726487
In [99]:
df_result_limits, df_result_xch, df_result_gamma = almacenar_resultados_modelo(prob)
df_result_limits['hora'] = pd.to_datetime(df_result_limits['hora'])
df_energia_demanda_limite = pd.merge(df_result_limits, D_h, on='hora', how='left')
df_energia_demanda_limite = df_energia_demanda_limite.astype({'hora': 'string'})

En algún momento del día la producción sube al máximo posible. ¿Qué momento es ese y por qué?

In [100]:
df_result_limits, df_result_xch, df_result_gamma = almacenar_resultados_modelo(prob)
df_result_limits['hora'] = pd.to_datetime(df_result_limits['hora'])
df_energia_demanda_limite = pd.merge(df_result_limits, D_h, on='hora', how='left')
df_energia_demanda_limite = df_energia_demanda_limite.astype({'hora': 'string'})
In [101]:
produccion_total_limite = df_energia_demanda_limite.groupby(df_energia_demanda_limite['hora'])[['energía', 'demanda']].sum().reset_index()
produccion_total_limite['demanda'] = produccion_total_limite['demanda']/10

plt.figure(figsize=(18,8))

plt.plot(produccion_total_limite['hora'], produccion_total_limite['energía'], label = 'Producción de energía total por hora', marker='x')

plt.xlabel('Hora')
plt.ylabel('MWh')
plt.xticks(rotation=90)

plt.legend()
Out[101]:
<matplotlib.legend.Legend at 0x213fd514a40>
No description has been provided for this image
In [102]:
produccion_minima = produccion_total_limite['energía'].max()
hora_minima = produccion_total_limite['hora'].loc[produccion_total_limite['energía']==produccion_minima].values[0]
print(f'El momento en que la producción llega a su máximo posible es: {hora_minima}')
El momento en que la producción llega a su máximo posible es: 2023-07-01 20:00:00
In [103]:
demanda_minima = produccion_total_limite['demanda'].max()
hora_minima = produccion_total_limite['hora'].loc[produccion_total_limite['demanda']==demanda_minima].values[0]
print(f'El momento en que la demanda es mayor es: {hora_minima}')
El momento en que la demanda es mayor es: 2023-07-01 20:00:00

El momento del día en que la producción cae a su máximo posible es a las 8 PM, debido a que este es el momento con mayor demanda energética.

El coste es mayor que con la demanda predicha. ¿Cuál sería ese coste máximo?

In [104]:
status = prob.solve()
print("Status: ", LpStatus[status])
costo_limite = value(prob.objective) 
value(prob.objective)
Status:  Optimal
Out[104]:
7755835.807726487
In [105]:
dif =round((((costo_limite*100)/costo_inicial)),2)-100
print(dif)
72.02000000000001
In [106]:
print(f'El costo máximo para este escenario es {costo_limite}, comparado con el costo inicial de {costo_inicial_onoff}, lo que indica una subida del {dif}% en el costo')
El costo máximo para este escenario es 7755835.807726487, comparado con el costo inicial de 4497165.960429066, lo que indica una subida del 72.02000000000001% en el costo

Si sube todavía un poco más, el problema se convierte en no factible.

In [107]:
limite_viable = daily_demand * (coeficiente + paso)
prob, e_ch, x_ch, gamma_ch = hallar_limite_viable(limite_viable, C, H, P_c, e_ch, x_ch, O_c, S_ch, ON_c, gamma_ch, E_c, beta, alpha)
In [108]:
status = prob.solve()
print('Status:', LpStatus[status])
value(prob.objective)
Status: Infeasible
Out[108]:
7941831.5364006655

Apartado C (opcional)¶

Tanto en el apartado A como el B hemos optimizado la gestión de únicamente un día, pero tenemos datos de un mes completo. Si nos imaginamos que tenemos que planificar para todo el mes de julio completo, hay 2 estrategias posibles:

  1. Hacerlo a las bravas, es decir, formular el problema completo para todo el mes e intentar resolverlo.
  2. Hacerlo día a día, es decir, optimizar primero para el día 1 sin más, a partir de cómo se queda la situación al final del día 1 optimizar para el día 2 sin mirar más allá, a partir de ese resultado optimizar para el día 3, etc.

La ventaja de la primera aproximación es que obtendría la solución realmente óptima, pero tiene el inconveniente de que según va creciendo de tamaño el problema a resolver, el tiempo y recursos necesarios para resolverlo pueden terminar siendo excesivos.

En cambio, la ventaja de la segunda aproximación es que resolver el problema para cada día es rápido (como ya hemos visto), pero como cada día sólo optimizamos sin mirar más allá, tomaremos decisiones subóptimas en cuanto a coste conjunto, dado que es como si estuviéramos a ciegas de cuál va a ser la demanda del siguiente día en adelante.

Plantea y programa el problema para el mes completo (ya sea en la modalidad 1 o la 2) y comenta los resultados obtenidos.

En el caso de que te atrevas con ambas modalidades, compara los costes para ver hasta qué punto la modalidad 2 es subóptima con respecto a la 1.

Criterios de evaluación¶

Para las formulaciones matemáticas:

  • Corrección de las fórmulas
  • Cumplimiento con todas las especificaciones del problema
  • Nomenclatura intuitiva de las variables
  • Explicación de qué representa cada una

Para la traducción a código de las formulaciones:

  • Corrección del código y ejecución sin errores
  • Resultados acorde con lo esperado
  • Dar respuesta a las preguntas formuladas (aquí se pueden añadir al notebook tantas celdas de texto como se consideren necesarias, además de las del código propiamente dicho)
  • Buen estilo de código (estructuración, modularidad, documentación, etc.)

Para la parte opcional:

  • Puntos extra acorde con lo realizado